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Abstract 

The Kuramoto phase diffusion equation is a nonlinear partial dif- 
ferential equation which describes the spatio-temporal evolution of 
a phase variable in an oscillatory reaction diffusion system. Synchro- 
nization manifests itself in a stationary phase gradient where all phases 
throughout a system evolve with the same velocity, the synchroniza- 
tion frequency. The formation of concentric waves can be explained 
by local impurities of higher frequency which can entrain their sur- 
roundings. Concentric waves in synchronization also occur in hetero- 
geneous systems, where the local frequencies are distributed randomly. 
We present a perturbation analysis of the synchronization frequency 
where the perturbation is given by the heterogeneity of natural fre- 
quencies in the system. The nonlinear ity in form of dispersion, leads 
to an overall acceleration of the oscillation for which the expected 
value can be calculated from the second order perturbation terms. 
We apply the theory to simple topologies, like a line or the sphere, 
and deduce the dependence of the synchronization frequency on the 
size and the dimension of the oscillatory medium. We show that our 
theory can be extended to include rotating waves in a medium with 
periodic boundary conditions. By changing a system parameter the 
synchronized state may become quasi degenerate. We demonstrate 
how perturbation theory fails at such a critical point. 



I. Introduction 



The formation of spatio-temporal patterns is ubiquitous in natural and ar- 
tificial complex dynamical systems |T1 [21 [3]. In oscillatory media pattern 
formation is tightly connected to the process of synchronization and plays 
an important role in a variety of systems far from equilibrium, such as ar- 
rays of Josephson junctions |4j, the Beluzov-Zhabotinsky reaction [5],|6], car- 
diac tissue [7], neural systems [8] and spatially extended ecological systems 
O [ini [H] ■ Different mechanisms for pattern formation are known. One of 
these is the interplay between attractive interaction, e.g. diffusion, which 
mediates long range correlations, and heterogeneity, or disorder, driving the 
system away from a uniform state. However, while large amounts of spatial 
heterogeneity describe the reality of most natural and biological systems, not 
much about the pattern formation and synchronization in disordered oscilla- 
tory media is known. 

Synchronization in the sense of a mutual adjustment of internal frequen- 
cies |il2j does not necessarily imply a total reduction of the system dimension 
to that of a single component, i.e. completely uniform dynamics. Instead, 
even in the synchronized state parameters like phase can vary across the sys- 
tem while the phase differences remain bounded or locked. In that case one 
can define waves that travel along a phase gradient [13]. If the wavelength 
is smaller than the diameter of the system these waves are perceived as time 
periodic spatial patterns. Such waves are a prominent feature in regular 
low-dimensional reactor topologies of chemical oscillating reaction-diffusion 
systems [H]. Wave propagation in oscillatory systems, although experimen- 
tally more difficult to assess, is also observed and of much relevance in a 
biological, medical, ecological and epidemiological context [Tf [Q], \T0\ 

Beside spiral waves and turbulence, concentric ring waves patterns are one 
of the most prominent features in two dimensional oscillatory media. They 
are usually associated with the presence of local impurities in the system 
[HIS]. These pacemakers change the local oscillation frequency and are able 
to entrain their surroundings, which finally results in regular ring or tar- 
get patterns. However, concentric waves of surprising regularity occur also 
in heterogeneous systems, where the local frequencies are distributed ran- 
domly. This was first reported and explained in [12] and subsequently also 
observed in chaotic phase coherent systems [U [TT] . In [TB] it was shown 
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that a phenomenological description can be obtained with phase equations 
using the Kuramoto model fT]. The analysis revealed that the random na- 
ture of the medium itself plays a key role in the formation of the patterns. 
In order to utilize or control these patterns, in general, it will not be suffi- 
cient to understand the mechanisms leading to pattern formation. Of equally 
importance is a knowledge about the time and length scales involved [T7]. 
However, in disordered systems no analytic formula for such quantities, let 
alone the full phase profile, are known. Here we show how estimations can 
be obtained by perturbation theory. 

The objective of this paper is to derive first and second order perturbation 
terms for the synchronization frequency in the nonlinear Kuramoto Phase 
Diffusion equation (KPDE) given a time-independent distribution of frequen- 
cies in the system. We present two approaches to this problem. The first 
approach, as described in Section II, is based on a direct perturbation expan- 
sion of the KDPE. It includes as a special case, for vanishing nonlinearity, the 
exactly solvable situation of the inhomogeneous heat equation. The second 
approach presented in Section HI is based on a Cole-Hopf transformation of 
the KPDE to a stationary Schrodinger equation for a particle in a disordered 
potential. Here, classical Schrodinger perturbation theory can be applied and 
leads to the same expressions as the first approach. We extend the second 
approach to describe variations in systems with topological charges. Such 
solutions exist for system topologies with periodic boundary conditions. In 
Section IV we apply the theory to simple topologies and calculate the first 
and second order perturbation terms of the synchronization frequency at the 
example of a d-dimensional medium with topological charges and also for 
the two dimensional surface of a sphere. Throughout we confirm our ana- 
lytic results by direct numerical simulations. We deduce the dependence of 
the synchronization frequency on the size and the dimension of the oscillatory 
medium and demonstrate that the second order perturbation term changes 
its scaling behavior at the critical dimension d = 2. Below that dimension it 
diverges with the system size and for d > 2 it diverges for small frequency 
correlation lengths. Finally in Section V, we take a look at regimes which 
can not be described by perturbation theory. In particular, we observe a 
discontinuous change of the location of a dominant pacemaker center. 

Let us start by reviewing the nonlinear phase diffusion equations [1]. A 
heterogeneous oscillatory reaction diffusion system may be described by its 
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full dynamics 

X = F(X,r) + V2DX(r) (1) 

where F(X, r) describes an oscillatory nonequilibrium reaction at a position 
r given the vector of reactant concentrations X and the diffusion V^DX(r) 
in the system, where D is a diagonal matrix of diffusion coefficients and the 
second order spatial derivative has to be applied component wise. Here 
we always assume that the local dynamics at the different locations are stable 
limit cycle oscillations. If the diffusive coupling only leads to small deviations 
from these limit cycles and if the medium is locally isotropic the system can 
be reduced [2] to the dynamics of phases ■(9(r) of the form 

^9(r) = cj(r) + V^^9(r) +7(V?9(r))^ . (2) 

These simplified phase equations define the nonlinear Kuramoto Phase Diffu- 
sion equations (KPDE). They were introduced by Kuramoto in 1976 [18j and 
are obtained by the perturbative method of phase reduction, using averaging 
techniques, described in his seminal monograph from 1984 [1]. Here, uj{r) is 
the local natural frequency of oscillation, we have used a scaling of time to 
make the diffusion coefficient in front of the Laplacian differential operator 
equal to one and the parameter 7 controls the nonlinearity, or dispersion. It 
can directly be interpreted as the nonisochronicity, which is the shear rate of 
the phase flow near the limit cycle and describes the sensitivity of the phase 
velocities to changes in the oscillation amplitude [I]. 

The heterogeneity in the system may be parameterized by the sample vari- 
ance 

or some norm of the two point correlation function C(r, r') if the frequencies 
are random (but quenched), e.g. 

E[uj{r)uj{r')]-E[io^] =a^C{r,r') with ||C|| = 1. (4) 

For the phase equations (j2]) to be applicable to the problem Eq. ([1]) the 
relaxation time of the amplitudes must be small compared to the time scale 
of the phase evolution P], i.e. 177 -C 1. In the following we define u = arj 
with normalized frequencies rj and use cr as a parameter of the system. 
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For a simulation of the KPDE ([2]) on a discretization of the medium it is 
of advantage to use the discrete Kuramoto model pQ 



= u;„ + ^ Anm sin(?9m - + Bnm 7(1- cos(^9„ - ^n)) (5) 

m 

where the Laplacian of the medium is defined through the values Anm and 
the square absolute value of the gradient through the choices of -Bnm- On a 
square lattice with nearest neighbor coupling of spacing h we have Anm = 
R — 



II. Perturbation Approach 1 

In synchronization the phase velocities of Eq. ([2]) have adapted to a common 
frequency 

t9(r) = n = ar/(r) + V¥(r) + 7 [Vt9(r)]^ (6) 

and with ^V?? = W'd = the phase gradient becomes stationary. In a 
homogeneous system, without disorder o" = 0, the constant phase profile 
'i9(r) = '(9° = solves the KPDE in synchronization, Eq. (jS]), with i7 = = 0. 
In contrast, in the presence of heterogeneity a > it is hard to obtain the 
stationary phase profile because the synchronization frequency is not known 
and must be calculated self-consistently [16]. Here we follow a perturbation 
approach by expanding in powers of the disorder a. Thereby, as will be 
shown below, non-trivial results are obtained in the second order. 

Given the normalized frequencies ri{r) it is possible to derive the pertur- 
bation series 

n{a) = an^^^ + (T2fi(2) + 0(^3) (7) 
directly by inserting the ansatz 

^ = + + 0{a^) (8) 

into Eq. ([6]) and regrouping terms according to powers of a 

00 

Q = J2(^^ (L^?^-') + b(^)) . (9) 
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Here, L = is the Hermitian, Laplacian operator, i!)^^^ is the perturbation 
term of order j in Eq. and the functions b'^-'^ are given by 

5(1) (r) = r){r) 

6(2) (r) = 7[V^9«(r)]' (10) 
i-i 

5(i>i)(r) = 7 ^ V??^^ V'i?^^-*) . 

i=l 

To proceed, it is convenient to expand into eigenfunctions of the Laplacian 
L. For a medium of finite volume and appropriate boundary conditions the 
eigenvalues of L are discrete. With the orthonormal eigenfunctions 
of the Laplacian and using the inner product of complex functions f and g 
defined for all positions r G M of the medium 

(ft.g) = / drr{r)g{r) (11) 

we can define the projectors 

Po = Pop|, and Qo = I-Po (12) 

with the identity operator I and the constant function po{r) = 1/ a/|M| which 
is the normalized eigenfunction of L to the eigenvalue Eq = 0. The operator 
Qo removes, in fact, the average from a function. Applying these operators 
to Eq. we obtain 

= Pob(^) (13) 

= L'i?^^) + Qob(^') . (14) 
The inverse operator of L in the image space of Qo is 

We can thus solve Eq. (HM and find the perturbation terms t!}^^^ up to a 
constant phase shift as 



^o) = -j;l^— ^p, . (16) 
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The equations (fTUj [T^ [TB]) can be iterated to obtain the full perturbation 
series ([3, [HD up to arbitrary order. Using the identities 



(plV^-Vpfc.) = ~Ek6kk' (17) 

(pI ■ v) (v^ ■ Pfc') = mvk' (18) 

and after some algebra we find for the first order and the second order per- 
turbation term of the synchronization frequency 

= (^>S.ste. (19) 

The coefficients rjl are the square amplitudes of the k^^ spatial Fourier modes 
of the frequencies, with respect to the system Laplacian. For k these 
values do not depend on the mean value of r/(r). Note that for isochronous 
oscillations 7 = the terms i)^^^^^ = are zero and the phase profile in 
synchronization is given exactly hy ^ = a^'-^^ and Eqs. ([IOl[T6D. In that case, 
the phase diffusion equation ([2]) is linear and readily solved in the Fourier- 
Space. 



III. Perturbation Approach 2 

In this section we will re-derive Eqs. ( IT9|20l) from a different point of view 
and in a somewhat more general form. It is well known that a nonlinear 
Cole-Hopf transformation 

?9(r) = ilnp(r) (21) 

7 

changes the KPDE ([6]) into a linear equation 

n-f p{r) = [7(T r]{r) + V^] p{r) = - H p{r) (22) 

for the ground state po{r) = p(r) of a Hamiltonian H with diagonal dis- 
order, given by the frequencies —'yar]{r), and ground state energy —'jQ 
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(see, e.g. Considering the frequencies ?7(r) a perturbation of strength 

e = 70", Schrodinger perturbation theory will give exactly the same re- 
sults for the synchronization frequency as obtained in the previous section 
Eqs. fll9|20p . However, as will be shown below the stationary Schrodinger 
equation is only one from a family of linear problems which are equivalent 
to the KPDE. 

To understand this, notice that there are two pitfalls to the transforma- 
tion Eq. (12T|) . First, it is only defined for non vanishing values 7 7^ 0. And 
secondly, since the ground state po{r) = 1/ of the unperturbed system 
Eq. (12^ is unique one is tempted to believe that the same is true for the 
homogeneous phase profile ("i? = const) of identical oscillators in synchro- 
nization. However, this is not necessarily the case because the phases '(?(r) 
are elements of a circle while p(r) is a real number. If the phase changes 
along a closed path from zero to a multiple of 27c it is a continuous func- 
tion on this curve while p is necessarily discontinuous. Indeed, for periodic 
boundary conditions multiple stable synchronized states can exist [19] (we 
will give examples for such inhomogeneous solutions in the next Section). 

We will therefore not take 'i9°(r) = as in the previous section, but in- 
stead assume a general phase profile 'i9°(r) in stable synchronization. Note 
that it is always possible to divide the phases formally into a time indepen- 
dent 'gauge' field ^9° and a time dependent deviation (f from that gauge field, 
^{r,t) = '(9°(r) + (p{r,t) . This corresponds to local rotations (i.e., the gauge 
field i)^{r) defines a position dependent choice of the coordinate frame) and 
yields new phases v^(r), so that ip = where in the old frame = We 
can define 

f]0(r) = V¥° + 7(V?9°)' . (23) 

Then the KPDE of the full heterogeneous system in synchrony take the new 
form 

n = uj{r) + fi°(r) + VV + 27VWy5 + 7 (V(/?)^ (24) 
and the gauge modified Laplacian reads 

L = 27Vi9°V + . (25) 

After the Cole-Hopf transformation ip{r) = ilogp(r) we find the eigenvalue 
problem 

7^ p = [-f(Tr]{r) + fi°(r) + 27V?9°V + V^] p . (26) 
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The change of gauge is nothing else but a similarity transformation of the 
Hamiltonian and does not effect its eigenvalues, unless the gauge field itself 
contains topological charges. If we choose a synchronized solution of Eq. ([6]) 
as gauge field so that n°(r) = is constant Eq. fl26l) can be written as 

7(r]-fi°) p = [7aV^ + 27V#V + V^] p . (27) 

where a non-Hermitean operator L = 27V'i9°V + is perturbed by a di- 
agonal disorder eV^ = —e diag(r7) with strength e = ■ya. Note that the 
constant function po{r) = 1/ a/|M| is an eigenfunction of L to the zero eigen- 
value -Eq = 0. This corresponds to a constant phase shift from ■(?°(r), i.e. the 
stable synchronization manifold for which we seek a perturbation expansion. 
Given orthonormal left and right eigenfunctions and p^ of L and the 
corresponding eigenvalues E^, the terms in the perturbation series for the 
eigenvalue Eq of — H with the largest real part 

7 (fi - = Eo = eEi'^ + e^^f + 0{e^) (28) 

are found to be 



E«' = (Pj-V,p„ 



4'' = -E 



(29) 



The discrimination between left and right eigenfunctions is necessary be- 
cause L is not Hermitian unless V"*?*^ = 0. For certain regular topologies 
and symmetric solutions i}^ the left and right eigenfunctions P^ = p^ are 
identical, nevertheless. In this case, we obtain the perturbation terms for the 
synchronization frequency as 

Q = Q"" + aQ^^^ + aV> +0{a^) (30) 



= (^>S.ste. (31) 

= _^J_y4 (32) 
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with the spatial Fourier modes 7]^ of the frequencies. This result can directly 
be compared to Eqs. (17lll9|20l) . The first noticeable difference is, that the un- 
perturbed system may have a synchronization frequency fl^ which is different 
from zero. The second difference is more subtle. The eigenvalues may 
have non vanishing imaginary parts. This corresponds to oscillatory modes 
during the transient to synchronization. Nevertheless, the sum Eq. ( 13^ and 
perturbation expansion Eq. fl5U]) are real if eigenvalues and eigenfunctions 
occur in complex conjugated pairs. 



IV. Examples 

Solution in a rectangular medium 

In the following we are interested to apply these results to some simple 
topologies. In order to apply our perturbation approach the spectrum of 
the Laplacian has to be calculated for every topology of interest. We start 
by examining a simple lattice. Let us consider a d- dimensional oscillatory 
medium L'^ C M.'^ with periodic boundary conditions and quenched random 
frequency disorder. We first note, that a constant phase gradient = , 
where I = [li, . . . ,1^) is a d-dimensional integer vector of winding numbers, 
solves the homogeneous equation Eq. fl25]) with 

^' = i[yJ\1\' ■ (33) 

Because of the periodic boundary conditions of the medium and the phases 
= + 271 we can include topological charges without phase singularities. 
The phase gradient is bounded and the use of the Kuramoto phase diffusion 
equation is justified. The unperturbed operator, Eq. (!25|) . reads 

L = 27^/^V + V2 (34) 

and the left and right eigenfunctions and eigenvalues that fulfill the periodic 
boundary conditions coincide and are simple harmonics 

Pk(r) = L-ie'^'^'" (35) 
LPk = (z27it.k-|knpk = i^Opk . (36) 
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The vector k = {ki . . . , kd) is also an integer vector, labeling the Fourier 
modes in the various directions. Eq. (1321) for the second order correction of 
the synchronization frequency shift gives 

This sum over the (i- dimensional integer lattice is potentially divergent de- 
pending on the small wavelength behavior of the terms 77^. Delta correlated 
random frequencies lead to an ultraviolet divergence in dimensions d > 1 
larger than one. We therefore have to restrict the perturbation theory to 
cases, where nearby frequencies are correlated, for instance as 

d I - 

E[ry(r)r/(r')] -E[77]^ = {2'k\^)~^ e'^^ (38) 

and A is some correlation length. Note, that Eq. fl38l) can only be an ap- 
proximation for small correlation lengths compared with the system size L, 
disregarding boundary effects. The expected value of the Fourier coefficients 
for |k| 7^ is then 

E [r^i] = e-2(^)'|i'l^ . (39) 

Using this expression one can calculate the expected second order perturba- 
tion terms in Eq. (!37|l numerically. An exact analytical expression exists in 
the simplest case of (i = 1, / = and delta correlated frequencies with A = 0. 
Then expression ( |37I) becomes 

A:>0 

where we could use the property of the Riemann zeta-function C(2) = vr^/G. 
In Fig. [1] we compare the shift of the synchronization frequency due to het- 
erogeneity for one dimensional systems with periodic boundary conditions, 
different lengths and winding numbers /. In order to observe the second or- 
der terms the linear contribution to the frequency shift must be exactly zero. 
This is achieved by shifting the average frequencies to zero (^)gygtgj^ = for 
each realization. The second order perturbation term is not affected by this 
change into a co-rotating frame of reference. The figure confirms the asymp- 
totic behavior of the synchronization frequency for 7cr <^ 1 and we find the 
scaling relation 

VL-VLq^ 7(7^ (41) 
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Figure 1: Influence of frequency disorder in a one-dimensional lattice. Plotted is 
the shift of the synchronization frequency from the frequency of the synchronized 
state with identical oscillators, f^O) as a function of the oscillator heterogeneity a. 
We compare simulations of the discrete Kuramoto model Eq. (0) for a chain of 
N oscillators (crosses) with the second order perturbation theory (solid lines), 
Eq. ()37p with L ^ N. Each cross from the simulations is an average value from 
50 runs with different realizations of iid. random frequencies {Ckk = l), where for 
each realization the average frequency has been shifted to zero {'n)system ~ ^' -l^^ft' 
comparison of the results for a ring of = 100 nonisochronous (7 = 0.25) phase 
oscillators, without topological charge (blue graph) and with a topological charge 
of ^ = 10 (red graph). Right: comparison of rings of different sizes = 16 (red 
graph) and N = 100 (blue graph) oscillators but with the same nonisochronicity 
7 = 0.25 and topological charge / = 3. 

as was previously observed in [16]. 



The scaling of the second order perturbation term with the system size L 
and the correlation length A can be studied by approximating the sum with 
a d- dimensional integral over |k| > 1 



2-d 



dk 



|k|>l 



(42) 



7 L^-^ x^-' r 



d-2 



X 



with X = —'K\f2 
ij 



We have here omitted constant factors, for instance from the integration 
over the (i-dimensional sphere shells |k| = const or upper and lower bounds 
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that allow the estimation of the sum from an integral. Depending on the 
dimension d we can use different asymptotic scaling relations of the incom- 
plete gamma-function in Eq. P2l) for a; — > 0, i.e. large system sizes or small 
correlation lengths. We find 

E ~ O (L) for d=l 

E[f)(2)j _ o[\og(^^ for d = 2 (43) 

E[f)(2)] ~ ©(A^-'^) for rf>3 . 

The analysis for a rectangular medium with no-flux or open boundary con- 
ditions gives analogous results, with the only difference that no topological 
charges are possible, i.e. |/| =0. 



Solution on a sphere 

Of special interest may be the synchronization frequency of a heterogeneous, 
oscillatory, reaction diffusion system on the surface of a sphere as a model for 
catalytic surface reactions on spherical bodies. Unlike in the torus topology 
of a rectangular medium with periodic boundary conditions, on the sphere 
topological charges always occur in vortex pairs of opposite charge. The 
method of phase reduction is not applicable in the vicinity of such phase 
singularities which act as fast pacemakers for the system. We will therefore 
only study perturbations of the homogeneous synchronized solution = 0. 
The eigenfunctions of the Laplacian on a sphere of radius R are spherical 
harmonics with 

Pim = ^Yim with / = 0, 1, . . . and m = 
ri 

(44) 

If we assume a homogeneous, isotropic distribution of frequencies on the 
sphere, the frequency correlator must have an S'0(3) symmetry. Let Ur be a 
transformation with Urez = Cr, which first rotates a vector in z-direction e^, 
around the y-axis to the zenith of Br and subsequently around the 2;-axis to 
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(a) 



(b) 



Figure 2: Quasi-regular wave patterns (a) in a rectangular medium with periodic 
boundary conditions and (b) on the surface of a sphere. We have used the discrete 
Kuramoto model Eq. dS} on (a) a grid of 150 x 150 phase oscillators and (b) on 
an almost homogeneous discretization of the sphere surface with 20480 points on 
the faces of a triangular tessellation [2lj . The natural frequencies of the individual 
oscillators were independently uniformly distributed with standard deviation a = 
0.2. The topology of the square lattice in subfigure (a) is that of a 2-torus and 
it is possible to have topological charges without phase singularities (large phase 
differences). We used an initial condition with topological charges Ix = 'i and 
ly = 7. Shown is the sine of the phases in gray levels after a transient time to 
synchronization. 

its azimuth. If we make the ansatz for an isotropic, homogeneous correlation 
kernel 

1 ~ /o; I 1 
E [r,{v)r^{v')] - E [r,f = Q J F,o (U^-^r') (45) 

we find the expected values of the spherical harmonics spectrum of the 
quenched frequency disorder as 



If the frequencies //(r) are delta correlated all coefficients q with / > are 
equal to one and the sums 




(46) 




(47) 



in Eqs. ( l20f32l) are divergent. We, therefore, have to assume a cutoff at a 
wave number /max ~ R/\ where A is a correlation length. This cutoff can 
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be sharp or exponential, as in the previous example, and we obtain again a 
relation 

E[fi(^)]~o(^log(^^)) . (48) 



V. Failure of Perturbation Theory 

The approach to the synchronization problem based on Schrodinger pertur- 
bation theory has another, conceptual advantage. Upon the variation of a 
perturbation system parameter the eigenvalues of a Hamiltonian can become 
degenerate or quasi degenerate. The perturbation theory of finite order must 
fail before such a point. This effect is illustrated in Figure ([3]), which inves- 
tigates the difference equation 

Ek Pn = (TVnPn + Pn^l - '^Pn + Pn+1 , (49) 

an approximation of Eq. (!22l) with 7 = 1 on a one-dimensional lattice and 
with open boundary conditions. Two pacemaker regions of different size and 
natural frequency are competing as wave centers. For low natural frequencies 
the larger and slower pacemaker region is dominating. By increasing both 
frequencies by a common factor a, the smaller and faster region gains ad- 
vantage. Two centers of waves can coexist in a small neighborhood around 
a critical value acr- Since the ground state of a one dimensional Schrodinger 
equation cannot be degenerate for a potential with finite square integral norm 
the largest and the second largest eigenvalue never coincide. The levels can 
come exponentially close depending on the distance between the two poten- 
tial wells. While the location of the dominating wave center shifts quickly 
upon variation of a the transient time until dominance is established scales 
as \Eo — £"1!"^. In practice, since close to the critical parameter value only 
transient behavior can be observed, one cannot say how far the boundary of 
the concentric waves will shift in either direction. The same transition can 
occur in a heterogeneous system with random frequencies as illustrated in 
Figures (Figs. H|5l) . The question where the pacemaker region of a heteroge- 
neous oscillatory medium is located given the local frequencies a;(r) cannot 
easily be answered without fully solving equation Eq. fl26|) numerically. 
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VI. Discussion 



In our study we have investigated the nonhnear Kuramoto Phase-Diffusion 
Equation (KPDE) in synchronization. We have apphed perturbation theory 
to calculate the synchronized state in a heterogeneous oscillatory medium. To 
our knowledge we have presented the first explicit analytical results regarding 
the oscillation frequency and the phase profile in such a system. Further we 
have identified different scaling relations depending on the system size and 
dimension and the frequency correlation length. We have shown that the 
perturbation approach can straightforwardly be applied to simple topologies. 

The first two terms of the perturbation series Eqs. fl30ti32p are intuitively 
quite meaningful. If a medium with random frequencies synchronizes to a 
common synchronization frequency Q then one expects Q to be close to the 
mean frequency in the system. But this is exactly the first order perturbation 
term. Any deviation from the mean frequency is due to the nonlinearity 7 
which appears as a factor only for higher order perturbation terms. 

Solutions of the Schrodinger Equation in disordered media are known to 
exhibit localization transitions [20l [21] , depending on the system dimension 
and the strength of the disorder. Given equation Eq. fl22]) and the properties 
of the disordered potential, all the results from condensed matter physics 
[211 [22l [23] dealing with the localized states in the impurity band, and in 
particular the ground state, can in principle be applied. One of the results 
is that in one and two dimensions all states are localized. It is straightfor- 
ward to show that in the limit of infinite system size perturbation theory 
does, in fact, not work for localized states (cf. one dimensional delta poten- 
tial). However, a perturbation ansatz is justified for states with a localization 
length larger than the system size. We have shown in the examples that the 
perturbation terms scale and diverge with the system size in one and two 
dimensional media. The result Eq. P3l) suggests that d = 2 is the critical 
dimension for the scaling of the synchronization frequency with the system 
size and the frequency correlation length. In one and two dimensions our per- 
turbation theory only gives good quantitative predictions for finite systems 
with 7crL -C 1, where L is the length of the system. In higher dimensions 
the synchronization frequency exists in the thermodynamical limit of L — cxo 
but it scales with the correlation length A of the frequencies as A^"*^. 
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The two presented perturbation approaches are equivalent in the sense that 
they lead to the same expressions for the perturbation terms, but when ap- 
plied to a specific realization of frequencies it can be of advantage to choose 
one approach over the other. By reducing the nonlinear KPDE in synchro- 
nization to an eigenvalue problem one can find the ground state energy and 
the corresponding phase profile in the discretized system to arbitrary order 
precision using linear algebra methods. Special attention must be given, if 
the phase profile spans phase differences over several decades, i.e. when the 
system size is large compared to the wave length. Then the exponentially lo- 
calized ground state must be computed to high precision even in the regions 
where it is several hundred orders of magnitudes smaller than at the localiza- 
tion point. For small nonlinearity 7 this second approach has no advantage 
over using the perturbation method Eqs. (fTOj [T3| [T6|) on the KPDE directly. 
In particular, the nonlinear Cole-Hopf transformation fl2T|) introduces addi- 
tional numerical errors. 
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DFG through the SFB555 and the Volkswagen Foundation. 
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Figure 3: Synchronization analysis in a one dimensional chain of = 200 non- 
identical Kuramoto phase oscillators with two competing pacemakers regions ac- 
cording to Eq. (f49l) . (b) frequency values —rjn corresponding to the potential of 
the discrete Hamiltonian (values ry„ are shifted so that the mean ^77) 



System 







is exactly zero). There are two potential wells at which the ground state can be 
localized, a deeper well on the left and a shallower but broader well on the right, 
(a) and (c) largest eigenvalues (solid blue lines) of the negative Hamiltonian — H 
in dependence on the heterogeneity a, the second order perturbation approxima- 
tion (Eq. ([29l) . dashed red line), and the value acr for which the ground state 
becomes quasi degenerate (dashed black line). The quality of the approximation 
can be seen in double logarithmic scales in subfigure (c). (d)-(f) numerically de- 
termined groundstate eigenvectors in a semilogarithmic scale for (d) a = 0.03, (e) 
a = acr = 0.036917 near the point of quasi degeneracy and (f) a = 0.04. Expo- 
nential localization at a potential well corresponds to concentric waves around this 
pacemaker region in the Kuramoto phase diffusion equations. 
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Figure 4: Synchronization analysis in a one dimensional chain of = 200 non- 
identical Kuramoto phase oscillators according to Eq. (149^ and independent, iden- 
tically, uniformly distributed random frequencies rjn (values are shifted so that the 
mean {v) System — ^ exactly zero), (b) frequency values —rjn corresponding to 
the potential of the discrete Hamiltonian (solid blue line) and a Gaussian filter- 
ing of width 2 (bold red line). There are several potential regions at which the 
ground state could be localized, (a) and (c) largest eigenvalues (solid blue lines) 
of the negative Hamiltonian — H in dependence on the heterogeneity a, the second 
order perturbation approximation (Eq. ([29l) . dashed red line), and the value acr 
for which the ground state becomes quasi degenerate (dashed black line). The 
quality of the approximation can be seen in double logarithmic scales in subfigure 
(c). (d)-(f) numerically determined groundstate eigenvectors in a semilogarithmic 
scale for (d) a = 0.04, (e) a = = 0.065135 near the point of quasi degener- 
acy and (f) a = 0.07. Exponential localization at a potential well corresponds to 
concentric waves around this pacemaker region in the Kuramoto phase diffusion 
equations [HI [15] . 
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Figure 5: Synchronization analysis in a one dimensional chain of = 200 non- 
identical Kuramoto phase oscillators with periodic boundary conditions and inde- 
pendent, identically distributed normal random frequencies rjn (values are shifted 
so that the mean (^)gyg^-£j^ = is exactly zero). The topological charge is / = 4. 
(a) three largest eigenvalue real parts ReEo > ReEi > ReEi of the difference 
operator in Eq. (j49p under variation of the heterogeneity a. The second and third 
eigenvalue in this example are complex conjugated for a < 0.09, and the ground 
state becomes quasi degenerate for ~ 0.154. (b) unperturbed {a = 0) rotat- 
ing wave solution on a ring of oscillators, here as the sine of the phase (in grey 
levels) on the side of a cylinder for illustration, (c, d) stationary phase profiles of 
the corresponding discrete KPE (Eq. ([5]), blue lines) for a = 0.15 and a = 0.16 
respectively, and the dashed black line is the unperturbed constant phase gradient. 
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